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We present linear response theories in the continuum capable of describing continuum spectra and dynamical 
correlations of finite systems with no spatial symmetry. Our formulation is essentially the same as the continuum 
random-phase approximation (RPA) but suitable for uniform grid representation in the three-dimensional (3D) 
Cartesian coordinate. Effects of the continuum are taken into account by solving equations iteratively with a 
retarded Green's function. The method is applied to photoabsorption spectra in small molecules (acetylene and 
ethylene) and inelastic electron scattering from a deformed nucleus 12 C. 
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CN ' Introduction 



Drip-line nuclei are weakly bound finite fermion systems. One 
would naturally expect that the continuum should be taken 
into account explicitly in the description of their excited 
states. Bound (L 2 ) representation, which is a popular pre- 
scription to construct spatially localized states, may give the 
gross properties of these states embedded in the continuum. 
However, the bound representation is incapable of describing 
quantities such as escape width and continuum spectra. The 
inclusion of the single-particle continuum for particle-hole (p- 
h) excitations has been achieved by using the so-called contin- 
uum random-phase approximation (RPA). 1 ' The continuum 
RPA combined with the Skyrme Hartree-Fock (HF) theory 
has been extensively utilized to study giant resonances (GRs) 
in spherical nuclei 2 ' 3 - 1 . Photoabsorption in rare gases was also 
studied by using the time-dependent local density approx- 
imation (LDA). 4 ' 5 ' One drawback of the method is that its 
applicability is limited only to systems with spherical symme- 
try. In this paper, we would like to present a prescription to 
treat the full three-dimensional (3D) continuum in the RPA 
level. 

For finite electron systems such as molecules, theoretical anal- 
ysis of photoabsorption spectra above the first ionization 
threshold requires continuum electronic wave function in a 
non-spherical multicenter potential. Advances in measure- 
ments with synchrotron radiation and high-resolution electron 
energy loss spectroscopy have enabled us to obtain oscillator 
strength distribution of an entire spectral region originated 
from valence electrons. 6 ' The photon energy dependences of 
molecular photoelectron spectra have also been measured. 
The spreading width is less important in molecules than nu- 
clear GR cases; however, there is a small width associated 
with molecular vibrations and dissociations. 

Since the Hamiltonian in the Skyrme HF theory for nuclei 
and the one in the Kohn-Sham (KS) theory of the LDA for 
electron systems are almost diagonal in coordinate represen- 
tation, a grid representation in the coordinate space provides 
an economical description. A uniform grid representation in 
the 3D Cartesian coordinate is becoming a standard method 
for calculation of the ground state. 7 ' 8 ' We would like to use 
the same 3D basis for calculation of the continuum response. 
The main problem lies in how to incorporate the scattering 



boundary condition in the uniform grid representation. 

Our method is based on the modified Sternheimer method. 9 ' 
It recasts the linear response problem into solving the static 
Schrodinger-type equation with a source term. The problem 
is then to solve the equation with an appropriate outgoing 
boundary condition. Our recipe here is to solve iteratively 
the equation taking into account the boundary condition em- 
ploying a Green's function of a spherical potential, separat- 
ing the self-consistent potential into long-range spherical and 
short-range non-spherical parts. This is reviewed in the next 
section. 



Linear response in the continuum of 3D real space 

Exact treatment of the continuum is possible with the use 
of a Green's function with an outgoing boundary condition. 
For spherical systems, the Green's function can easily be con- 
structed by making a multipole expansion and discretizing 
the radial coordinate. This is an essential part of the contin- 
uum RPA method. In this section we present a method to 
construct a Green's function in the 3D grid representation for 
a system without any spatial symmetry. 

The linear response theory is formulated most conveniently 
using a retarded density-density correlation function. 10 ' 

n(r,r';w) = f die*"' - "^^, t; r', 0), (1) 

zU(r,t;r',t') = 9{t - t'){0\\p(r, t),p(r', t')]\0). (2) 

The RPA transition density in response to an external field 
T4xt can be expressed as 



5p{v,u) 



d 3 rn RPA (r, r'; w)Kxt (r, u) 
dVn (r,r';w)V sc f(r',u;), 



(3) 

(4) 



where the Ho is the independent-particle density-density cor- 
relation function which is defined by identifying the state |0) 



in Eq. (0) with the unperturbed HF (KS) ground state and 
assuming that the density operator is evolved in time with the 
static HF (KS) Hamiltonian. V sc f, the self-consistent field, is 
the sum of the external field and the dynamical induced field: 



where functions ipi are defined by 



i>i(r;E,V scl )= / dVG i+ V,r';£)V; cf (r')<Mr). (11) 
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Using properties of the Green's function, we have an equation 



where V[p(r)] is a single-particle potential of either the HF or 
the KS equation. In order to treat the continuum, we express 
the no in the form 1 ' 4 -' 



no(r,r';u,) = ^ &(r) { (G (+) (r,r' ; (e, - u)*))* 

i 

+ G (+) (r,r';e 1+ ^)}^(r'), (6) 

assuming that the occupied states have real wave functions 
4>%. G' + ' is a single-particle Green's function with an outgoing 
asymptotic boundary condition, and is defined by 



G (+) (r,r';£) = (r\(^E+^V 2 -V(r)+ir^ ^r'p) 



E-e h +m' 



(e + ^-V 2 - V\ M^ E, V acl ) = K* (r)cfc(r). 



(12) 



(8) 



The integral in Eq. (nil) is thus converted into a differen- 
tial equation Q13) . This procedure is known as the modified 
Sternheimer method. 9 ' To determine a solution of Eq. (113) 
uniquely, we have to specify the outgoing boundary condition. 
This is our final task. 

We assume that the potential V can be divided into two 
parts, a long-range spherical part Vo(r) and a short-range 
non-spherical part v = V — Vo(r). For instance, Vo(r) is taken 
as a Coulomb potential induced by a spherical charge distri- 
bution. A necessary condition is that v must vanish outside 
the box (model space). Let us denote a Green's function for 
the Hamiltonian Ho = -V 2 /2m + Vo as G^(E), and first 
calculate an action of G on a vector $, ^(E) — Go(E)&, 
by solving the differential equation 



where <j>k and tk are the HF (KS) single-particle orbitals and 
their eigenenergies, respectively. 



{e - (-^V 2 + Vb(r)) } *(r; E) = $(r). 



(13) 



In writing the Green's function in terms of single-particle or- 
bitals, as in Eq. (M), we need to calculate both occupied and 
unoccupied orbitals including the continuum wave functions 
with positive energies. This would be a difficult task. The 
trick of the continuum RPA is an explicit construction of the 
Green's function in the coordinate space. If the potential V(r) 
has spherical symmetry, G^ + ' (r, r'; E) can be reduced into the 
radial function g\ \ (r, r';E). Then, we can easily construct 
9\ +) as 



The solution at the boundary of the box (r > R) is prepared 
using a multipole expansion method. 



*(r; E)1>r = ]T Wl { *"> E) Y lm (m m (E), 



(14) 



*, m (£) = 2m f d 3 r' Ul{r '; E) Y lm (r')f(r'), (15) 
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(9) 



where ui and w\ are independent solutions of the radial 
HF (KS) equation. The asymptotic boundary condition of 
g\ is given by an asymptotic behavior of the irregular solu- 
tion w\ (r). Unfortunately, this explicit construction of the 
Green's function is possible only for a one-dimensional case. 
It is difficult to impose the boundary condition for a case of 
a deformed potential. Furthermore, for 3D cases, the num- 
ber of spatial grid points is too large to construct the Green's 
function (iV gr id x iV gr id matrix). 

Instead of explicit calculation of the Green's function, we shall 
utilize a modified Sternheimer method. The transition den- 
sity can be written in the form 



where ui(r;E) and w ; (r; E) are solutions of the radial 
differential equation being normalized as the Wronskian 
W[w;,uij ] is unity. Then, we use the following identity for 
the Green's function 



G (+) (E) = G +) (E) + G^ } {E)vG m (E), 



(16) 



to calculate actions of G^ + '(E). Namely, instead of solving 
Eq. (jld) , we solve the equation 



(1 - G { +) (E)v)ME) = G +) (£)U scf ^, 



(17) 



by using an iterative method. In this way, we fix an outgo- 
ing asymptotic boundary condition for the Green's function 

G (+) (£). 



6p(u) = J2& {(<M(^ - w)*))* + ipi(u +uj)}, (10) 



Now let us describe a numerical algorithm for calculation of 
transition density Sp. There are three-nested iterative loops 



(I D II D III) to solve the multilinear equations: 

I) Solve Eq. (§); (1 - Ii ^)8p = IL.Kxt. 

II) Calculate actions of G^ by solving Eq. (M/fl) . 

III-l) Evaluate boundary values of actions of Gq using 
Eqs. (Q and (|l|). 

III-2) Calculate actions of G { +) by solving Eq. ([jj|). 
We use the generalized conjugate gradient method for non- 
hermitian problems, I) and II). For hermitian problems such 
as III-2), the conjugate gradient method is the most efficient 
algorithm. In the procedure I), if we neglect the dynami- 
cally induced field 5V/8p, the calculated response is called 
"independent-particle approximation" (IPA) in the following 
sections. 



Continuum RPA vs discrete RPA 

In this section, we compare results of the continuum RPA cal- 
culation with those of the RPA discretized in a spherical box. 
The modified Skyrme force, so-called "BKN interaction", 11 ' 
is adopted. 
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/iHF = 



^V 2 + - 4 t p(r) + -t 3 p(rf 



4 / r — r 



TT^). 



(18) 



where spin-isospin degeneracy (each nucleon with a charge 
e/2) is assumed. Thus, p(r) = 4^\ \<fii(r)\ 2 . The same pa- 
rameters as those in Ref. 11) are used in the calculation. The 
self-consistent HF potential is calculated in the spherical box 
of radius it = 8 fin for the 16 nucleus. The single-particle 
energies for 0s and Op orbitals are —28.4 and —16.7 MeV, 
respectively. 



The strength function is defined by 
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Fig. 1. Calculated monopole strength distributions for 16 0. The red 
line indicates results of the continuum RPA while the blue (green) 
line is the "discrete" RPA with a smoothing parameter r=l (5) 
MeV. 



matrix of 80 x 80. This size of matrices can be easily handled 
and the matrix inversion in Eq. (E2j) can be done. How- 
ever, for 3D calculations in the following sections, since the 
size of the matrix becomes of order 10 000, we shall use the 
Sternheimer method with the outgoing boundary condition 
described in the previous section. The strength function is 
calculated with Eq. (EOT. 

In order to see the effects of the continuum, we calculate the 
"discrete" RPA response as well. This can be done by simply 



using a standing- wave solution w\ (r) (w\ (R) = 0) in Eq. 



,(o), 



(+)/ 



(Bj) instead of the outgoing solution w t (r). Alternatively, 
we may solve the modified Sternheimer equation (113) with a 
boundary condition ipi = for r > R. 



S{E) = ^|(n|K x t|0)| 2 5(£-£ n ), 

= -hmJdrVl t (r)Sp(r;E), 
= -ilmTr^nRPAO^Kxt) 



(19) 

(20) 
(21) 



Since the obtained HF potential of O is spherical, we can 
calculate the continuum response in the radial coordinate fol- 
lowing the conventional continuum RPA procedure. 1 ' The 
single-particle Green's function is calculated as Eq. (BJ) and 
the RPA response function can be explicitly constructed by 



Calculated strength functions for isoscalar giant monopole 
resonance (ISGMR) are shown in Fig. hi The external field 
is Vext = r 2 in this calculation. For discrete RPA calculation, 
since all excited states are discrete, we use complex energies 
E — E — iT to produce a continuous profile of the strength 
functions. As we see in the figure, the continuum response 
is significantly different from the smoothed discrete response. 
In the discrete calculation, there are peaks characterized by 
the size of the box. In Fig. pi, concentrated strengths appear- 
ing around 40 MeV are spurious peaks produced by the finite 
box. The positions of these spurious peaks will change as the 
size of the box is modified. In the continuum calculation, a 
dip emerges at the threshold energy of the 0s orbital (28.4 
eV) but there is only a smooth tail beyond 30 MeV. 



n RPA (£) = 



5V 

Sp 



n (E). 



(22) 



The strength function is calculated with Eq. (E1J) . The 
adopted box is a sphere of radius it = 8 fin discretized in 
meshes of 0.1 fm. Thus, the response function n(r, r'; E) is a 



Inelastic electron scattering from C 

In this section, we discuss a 3D calculation of the continuum 
response. As we have mentioned above, although the contin- 
uum RPA has an ability to treat the single-particle continuum 
exactly, it has never been applied to deformed nuclei because 
of the difficulty of constructing the Green's function. There- 



fore, the continuum effect on deformed nuclei has rarely been 
investigated so far. We use the same BKN interaction as in 
the previous section and try the first calculation for a light 
nucleus 12 C. The model space is a 3D spherical box of radius 
R — 8 f m with square meshes of Ax — Ay = Az — 1 fm. The 
obtained HF ground state has an oblate shape in which the 
axis ratio is approximately three to two. 

Utilizing the Born approximation, the cross section of the 
inelastic electron scattering for momentum transfer q and en- 
ergy loss E is approximated by a product of the strength 
function and the Mott cross section. 10 ' 
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S(E,q), 



(23) 



where S(E, q) is defined by Eqs. (|19|-J2l|) with V cxt = e lqr (l + 
T3)/2J 1 J Since the ground state of ^C is deformed, the re- 
sponse to K x t depends on the direction of momentum transfer 
q relative to the orientation of the nucleus. We show results 
for two cases that the q is parallel and perpendicular to the 
symmetry axis (z-axis) of the nucleus. Of course, in experi- 
ments, the response averaged over different angles is observed. 

In Fig. H, the strengths as functions of energy are displayed 
for low momentum transfer (q — 0.4 fm -1 ) and for high mo- 
mentum transfer (q — 4 fm -1 ). It seems that peaks related 
to giant resonances can be seen up to 25 MeV. For the case of 
high momentum transfer, the response strongly depends on 
the direction of q relative to the intrinsic orientation. 



Photoabsorption of molecules 

First we discuss valence-shell photoabsorption of acetylene 
and ethylene molecules. The KS Hamiltonian for valence elec- 
trons is 



h KS = "^V 2 + tfon + e 2 fd 3 r'-0^ + Mxc [p(r)],(24) 



where Vi on is an electron-ion potential for which we em- 
ploy a norm-conserving pseudopotential 12 ' with a separable 
approximation. 13 -' The /j xc is an exchange-correlation poten- 
tial. 

It is well known that the energy of the highest occupied molec- 
ular orbital (HOMO) does not coincide with the first ion- 
ization potential in the simplest local-density approximation. 
This fact causes a serious problem for the continuum response 
calculation that the ionization threshold cannot be adequately 
described by the static KS Hamiltonian. Furthermore, the 
excited states around the ionization threshold appear in too 
low excitation energies. To remedy this defect, a gradient cor- 
rection for the exchange-correlation potential has been pro- 
posed. We utilize the one constructed by van Leeuwen and 
Baerends 14 ' which we denote as // -\ It is so constructed 
that the potential has a correct — e 2 /r tail asymptotically. 
The energy of the highest occupied orbital also approximately 

Because of the assumption of the BKN interaction, the 
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Fig. 2. Response of C to an electromagnetic probe. The red 
and blue lines correspond to momentum transfer q=0.4 and 4 
fm , respectively, (a) The momentum transfer q is parallel to 
the symmetry axis of the nucleus, (b) The q is perpendicular to 
the symmetry axis. 



coincides with the ionization potential. For small molecules, 
TDLDA calculations with this gradient correction have shown 
to give an accurate description of discrete excitations in small 
molecules. 15 ' 

In the following sections, we employ a sum of the exchange- 
correlation potentials of /v ' of Ref. 16) for the local-density 
part and ^jS lb ' for the gradient correction. We should remark 
here that an accurate calculation of the gradient correction 
// LB '(r) becomes difficult at far outside the molecule, because 
/i' LB ' depends on |Vn(r)|/n(r) 4 ' 3 which approaches a finite 
value but both the numerator and the denominator approach 
zero at r — -> oo. Thus, we use an explicit asymptotic form, 
// LB '(r) = — e 2 /r for r > R c . In the following, R c is chosen 



as 5.5 A. The spherical box is taken as radius R 
meshes of Ax = Ay = Az — 0.3 A. 



6 A with 



isoscalar response is shown in Fig. YA with V c 



The acetylene molecule, C2H2, has a symmetry configuration 
of Daah- This high symmetry has enabled calculation of the 
Green's function using a single-center expansion. 17 ' Even so, 
only two kinds of molecular orbitals, 3o g and ln u , which are 
primarily derived from atomic p states, have been taken into 
account in Ref. 17), because it was difficult to describe the 
s-derived states in the single-center formulation. Since our 
framework is free from this problem, we consider all valence 
orbitals, including the 2<r 9 and 2o u in addition to the above 
p-derived orbitals, to calculate the photoresponses. 



All atoms are located on the z-axis at ±0.601 A for carbon 
and ±1.663 A for hydrogen. There are ten valence electrons 
and the calculated energies of occupied orbitals are listed in 
Table 111. We obtain the photoabsorption oscillator strengths 
as a function of photon energy, as shown in Fig. H. The cal- 
culation indicates a sharp bound resonance at lu — 9.6 eV 
and a broad structure around 15 eV which seems to be a 
superposition of three resonances. The resonance at 9.6 eV 
strongly responds to a dipole field parallel to the molecular 
axis. The large oscillator strengths in the IPA at lu = 5 ~ 8 
eV and at lu — 12.5 ~ 13.5 eV are shifted to higher ener- 
gies by the dielectric correlation. The agreement with the 
experimental data is significantly improved by the inclusion 
of the dynamical induced field (screening). The static dipole 
polarizability is also affected significantly: In the IPA calcula- 
tion, the polarizabilities parallel (an) and perpendicular (aj_) 
to the molecular axis are ay = 10.7 A 3 and ai = 3.87 A 3 . 
The dynamical screening reduces these values to 4.79 A 3 and 
2.77 A 3 , respectively, which well agree with the experimental 
values, «n = 4.73 and a± = 2.87 A 3 . 18) 

Another example we show here is ethylene, C2H4. Ethy- 
lene is the simplest organic 7r-system, possessing the D2/1 
symmetry. Its two carbon atoms lie on the x-axis at 
x = ±0.6695 A and its four hydrogen atoms lie in the 
x-y plane at (x,y) = (1.2342,0.9288), (1.2342,-0.9288), 
(-1.2342, 0.9288), (-1.2342, -0.9288) in units of A. There are 
12 valence electrons, and calculated eigenenergies of occupied 
orbitals are listed in Table hi 

The calculated photoabsorption oscillator strengths are shown 
in Fig. m The agreement with an experiment 19,21 ' is excel- 
lent. Almost all the main features of photoabsorption spectra 
are reproduced in the calculation. The observed bound ex- 
cited states show different photoresponses according to the 
direction of the dipole field. The lowest peak at u — 7.6 eV 
is mainly a response to a dipole field parallel to the molecu- 
lar (C— C) axis. This is associated with the excitations of the 
HOMO lfc3 U electrons. On the other hand, states at 9.8 eV re- 
spond almost equally to a dipole field of x, y, and z directions, 
to which both 1633 and 163^ occupied orbitals contribute. A 
small peak at lu = 11.4 eV is calculated as a resonance of 
3a g orbital, which may correspond to a small shoulder in the 
experiment. Beyond 11.7 eV, the HOMO electrons are in the 
continuum. The first prominent peak at 12.4 eV is a reso- 
nance with respect to a dipole field of y direction which is 
in the molecular plane and perpendicular to the C— C bond. 
The excitations of lb2 U and lb3 g electrons are the main com- 
ponents of this resonance. In the region of 13.2 < lu < 20 
eV, 1&3 9 electrons can be excited into the continuum and pro- 
duce the smooth background of oscillator strengths (0.1 ~ 0.2 
eV -1 ). The peak structure at 14.6 eV is produced by the ex- 
citation of l&2u electrons. The peak at 16.4 eV is constructed 
by excitations from the 2bi u occupied orbitals, while the ex- 
periment indicates the resonance at 17.1 eV. 

In Figs, fcfl and H, we show about 75 % of the Thomas- Reiche- 
Kuhn (TRK) sum rule for valence electrons. In other words, 
one-fourth of the TRK sum rule value lies in an energy range 
of over 40 eV. 



Table 1. Calculated eigenvalues of occupied valence orbitals and 
experimental vertical ionization potential (IP) in units of eV. The 
experimental data are taken from Ref. 20) for acetylene and from 
Ref. 21) for ethylene. 
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Fig. 3. Calculated and experimental photoabsorption oscillator 

strengths of acetylene. The red line is the calculation compared 
with synchrotron radiation experiment (green line) 22 ) and high- 
resolution dipole (e,e') experiment (blue line). 23 ' The dashed line 
is the IPA calculation without dynamical screening. 



Conclusions 

We have developed a method based on the time-dependent 
density-functional theory of investigating responses in the 
continuum for systems with no spatial symmetry. The 
method allows us to treat the escaping width and the dy- 
namical correlation effects self-consistently. We have shown 
the ISGMR in the continuum for ie O and compared the re- 
sults of the continuum and discrete calculations. The con- 
tinuum response cannot be reproduced by simply smoothing 
the discrete strengths. Perhaps, this is because the contin- 
uum level density is not properly taken into account in the 
discrete calculation. Then, we have shown a calculation of 
electromagnetic response function for a nucleus 12 C and pho- 
toabsorption spectra for molecules, acetylene and ethylene. 
The main difficulty is the heavy computational task, espe- 
cially for calculations of response near the threshold energies. 
Using a single CPU of a Fujitsu VPP700E supercomputer at 
RIKEN, the numerical calculation of photoabsorption spectra 
of ethylene in Fig. takes about 30 minutes at a single photon 
energy lu. Since we need to calculate the response to dipole 
field of the x, y, and z directions, it takes about 1.5 hours to 
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Fig. 4. The same as Fig. 3 but for ethylene. The experimental data 
are taken from Refs. 19721). 



calculate a single energy point. However, we have found that 
the inclusion of the small imaginary part in energy, E + iF, 
facilitates convergence of iterative procedures in the calcula- 
tion, r also plays the role of lowering an energy resolution of 
the calculations. Thus, we can choose a value of T depending 
on the energy resolution required in each problem. 

The method is capable of calculating response functions of 
many-particle systems, below, near and above the separation 
energy (ionization threshold) in a unified manner. The results 
seem very promising and encourage us to apply the method 
to responses of drip-line nuclei in the near future. 
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